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o 

^ . One of the main interests in seismology is the formulation of models able to describe the clustering 

' in time occurrence of earthquakes. Analysis of the Southern California Catalog shows magnitude 

clustering in correspondence to temporal clustering. Here we propose a dynamical scaling hypothesis 
in which time is rescaled in terms of magnitude. This hypothesis is introduced in the context of 
a generalized trigger model and gives account for clustering in time and magnitude for earthquake 
occurrence. The model is able to generate a synthetic catalog reproducing magnitude and inter-even 
time distribution of thirty years California seismicity. 
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The great interest in the study of earthquake occurrence is linked to the challenge of predicting the time, the 
location and the energy of the next earthquake. The energy release E in a seismic event can be expressed by the 



o 

magnitude M via the logarithm relation M oc log E |1| , and the magnitude distribution is described by an exponential 
law usually referred as the Gutenberg- Richter (GR) law 0] P(M) ~ \Q~ bM , where b is a parameter close to one. 
. The logarithm relation leads to a power law behaviour for the energy distribution, which is generally the signature of 
' ^ ' critical phenomena. 

£>v It is widely observed that earthquakes tend to occur in bursts. These bursts start immediately following a large 
main event, giving rise to the main-aftershock sequences, described by the Omori law 0]. This states that the number 
of aftershocks n(t) decays in time as n(t) ~ (t + c)~ p where p is generally close to 1 and c is an initial time introduced 
in order to avoid the divergence at t = 0. The most important implication of this law is that we cannot assume a 
Poissonian occurrence for earthquakes, namely characterized by a constant rate of occurrence, but rather a clustered 



one. 

Another signature of non-Poissonian behaviour for earthquake occurrence is the complex distribution of the inter- 
occurrence times between two successive events. In fact, for a Poissonian process, this distribution would be an 
exponential whereas experimental data exhibit a more complex behaviour Q. Moreover, one can compute the intertime 
distribution D(At, Ml) where At is time distance between successive events occurred inside a finite geographic region 
and with magnitude greater than a given threshold Mr,. Indicating with Pc{M) the cumulative magnitude distribution 
inside the considered region, one observes 0, 

D(At, M L ) = P c (M L )f(P c (M L )At) (1) 

• i-H 

where / is a universal function, independent on Ml and on the geographical region. The observed universality is a 
further signature of criticality and indicates that D(At, Ml) is an appropriate quantity to characterize the temporal 
O r clustering of earthquakes. 

A widely used approach to earthquakes clustering is provided by "trigger models" [6j. These assume a Poissonian 
occurrence of triggering events, whereas the occurrence of the "triggered" earthquakes is described in terms of a 
correlation function with previous events. Among the trigger models the Epidemic Type Aftershocks Sequence 
(ETAS), introduced by Kagan-Knopoff and developed by Ogata describes mainshocks and aftershocks on the 
same footing. More precisely, each earthquake can generate "its own aftershocks" and furthermore the number of 
these aftershocks depends exponentially on the magnitude of the " main" . The model has been deeply investigated 
analytically and numerically [9j . 

In this paper we are interested in the description of temporal evolution of seismic activity. For this reason we 
neglect spatial dependencies and treat seismicity as a stochastic process Mj(tj), where Mi is the magnitude of the 
i — th earthquake occurred at time U inside a large but finite geographic region. The process is defined by the 
conditional probability density p(M(i)\{Mi(ti)}) to have an earthquake of magnitude M at time t given the history 
of past events {Mj(ij)}. Here we consider a generalized version of the trigger model by Vere- Jones 

p(M(t)\{Mi(ti)}) = p(M(t)\Mi(ti)) + nP(M) (2) 

i:ti<t 
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FIG. 1: (Color online) (a) Experimental distributions of the quantities Qj (red o) and Rj (black continuous line) for the 
California Catalog. Rj is measured in unit of (Ghours)" 1 and in order to improve the comparison Qj is vertically shifted by 
the constant amount —1.5. Peaks for Rj indicate main-after shock sequences, (b) Numerical distributions of Qj and Rj from 
Eq.0. 

where p(M (t)\Mi(ti)) is the "two-point" conditional probability density, /i is a Poissonian rate and the magnitude 
distribution P(M) ~ 10~ bM obeys the GR law. Different forms of p(M (t)\Mi(ti)) correspond to different models for 
seismicity. In the ETAS model one assumes || 



where the propagator g(t — tf, Mi) oc 10 aMi (t — ti + c)~ p . In order to have a normalized probability one must im pos e 
p > 1. Moreover, if a > b the model presents finite time singularity unless one assumes a large magnitude cut-off |lCj. 
Alternatively, one must take a < b as supported by some experimental observations 

A strong assumption of the ETAS model is the factorization in Eq.©, which states that the magnitude of an 
earthquake is completely independent on the magnitudes and times of occurrence of previous events. In order to 
test this assumption with real seismic data, we observe that the quantity Pc(Ml) ~ 10~ bML takes the role of a 
characteristic time scale in Eq.Q. Hence, if one considers a subset of TV events, the quantity Q = N/[^2 t 10~ fcMi ] 
can be related to the rate of occurrence R = — ij-i)], where the sum is inside the chosen subset. To this 

extent, we divide data recorded in the Southern California Catalog (1975-2004) |l2( in subsequent sets of 200 events 
with M > 2.5 and we compute the quantities Qj and Rj inside the j-th subset. If the magnitude distribution were 
constant in time, as supposed in Eq. Qj should fluctuate around an average value. Conversely, the experimental 
Qj displays scattered and narrow peaks (Fig. la). Interestingly, these peaks are closely located to peaks in the Rj 
distribution. It is well known that peaks of Rj are located soon after main-shocks and indicate the presence of 
main-after shock sequences. Fig. la, then, shows that in subsets of the catalog where activity has an higher rate, 
the probability to have large magnitude events is also raised. This aspect can be directly investigated by computing 
the cumulative magnitude distribution Pc(M) only inside the ensemble of main-aftershock sequences. Considering 
only sequences with main-shock magnitude M > 6, one obtains that Pc(M) exhibits a GR behaviour with a best 
fit 6-value b = 0.75, lower than the 6-value obtained for the whole catalog (b = 0.95) within the 95% significativity 
level. This result further supports the idea that large earthquakes not only produce the clustering in time described 
by the Omori law, but also a clustering in magnitude. The ETAS model does not take into account this last physical 
mechanism. 

In order to include the magnitude clustering within a trigger model approach, we propose a dynamical scaling 
hypothesis: the magnitude difference Mi — Mj fixes a characteristic time scale Tij = fclO^ J l> so that the 
conditional probability is magnitude independent when times are rescaled by Tjj and k is a constant measured in 
seconds 



p{M(t)\Mi(ti)) = P(M)g(t - U;Mi) 



(3) 



p(M i {t i )\M j {t j ))=F 



(U - tj) 



(4) 



Let us then consider the probability to have an event of magnitude M at time t given a triggering event at time to 
of arbitrary magnitude M , p(M,t- t ) = J dM P(M,t\M ,t )P(M ). Assuming the GR law for P{M ) and using 
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p(M, t-to) = 7—ttz / F{z^)dz. (5) 



Eq.Q, one finds 



(* - to) p Jo 

From this equation we obtain both the GR and Omori law independently of the specific form of F(z) provided that 
the appropriate constraints are imposed at small and large z. In fact, assuming that the conditional probability (6) is 
maximal soon after the triggering event, must be F(0) > 0. Furthermore, in order to have normalized distributions, 
the conditional probability must decay to zero for large time separation and a constraint on the behaviour of F(z 1 / p ) 
must be imposed at large z, namely a decay faster than 1/z. Because of this constraints, the integral in the rhs of 
Eq.JSJl is a constant for large t — to, and the GR and Omori law directly follows from Eq.QJ. The above observation 
suggests that statistical features of the trigger model can be independent on the detailed form of F(z) once the scaling 
Eq.Q is assumed. This hypothesis together with the relationship between numerical and experimental behaviour can 
be directly tested in numerical simulations. 

In a numerical protocol one assumes at initial time to — a single event of arbitrary magnitude chosen in a fixed 
range [M m i n , M max \. Time is then increased of a unit step t = to + 1, a trial magnitude is randomly chosen in the 
interval [M m j n , M max ] and Eq.(|2J) gives the probability to have an earthquake in the time window (to, to + 1). If this 
probability is larger than a random number between and 1, an earthquake takes place, its magnitude and time of 
occurrence are stored and successively used for the evaluation of probability for future events. Time is then increased 
and in this way one constructs a synthetic catalog of N e events. The term /x in Eq.@ represents an additional source 
of earthquakes Poissonian distributed in time with a magnitude chosen from the GR distribution with b = 0.8. 

Following this protocol, we generate sequences of 15000 events using a power law form for F(z) 

F{z) = - T ^— (6) 
■ J z x + 7 w 

and then we compute the numerical distributions D(At, Ml) and P(M). These distributions are compared with the 
experimental data from the Southern California Catalog. For different values of A, it is always possible to find a 
set of parameters A,j,b/p, \i such that numerical data reproduce, on average, the statistical features of earthquake 
occurrence both in time and in magnitude. The parameter k is fixed a posteriori in order to obtain the collapse 
between numerical and experimental data. 

In Fig. 2 we plot the experimental and numerical D(At, Ml) considering two different values of A (A = 1.2 and 5) 
and Ml (Ml = 1.5 and 2.5). In the inset we also present the magnitude distributions. Data for different values of 
the parameters follow a universal curve and the same collapse is obtained for other values of A > 1. The accordance 
between experimental and numerical curves indicates that the hypothesis of dynamical scaling is able to reproduce 
two fundamental properties of seismic occurrence, namely the GR law and Eq. Q), independently of the details of 
F(z) 0. 

The ETAS model is a particular case of Eq.© corresponding to 7 = and A ~ 1. We want to stress the important 
difference due to the presence of a non-zero 7. From a mathematicalpoint of view, the constant 7 avoids the finite 
time singularity of the ETAS model with a = b discussed previously |9| . From a physical point of view, the constant 
7 gives rise to the observed clustering in magnitude. Indeed, for a given mainshock of magnitude Mj at time tj, at 
each time (ti > tj) it is possible to define a sufficiently large magnitude difference AM such that, if Mj ~ Mi > AM, 
we have that z x is negligible with respect to 7 and therefore F[(ti —tj)/rij] ~ A/-f. In other words after a large event, 
small earthquakes tend to be equiprobable. 

We have also performed more extensive simulations using a different expression for F(z) 

e z — 1 + 7 

Eq.0 states that two events of magnitude Mj and Mj are correlated over a characteristic time r,j and become 
independent when f j — tj > Tij . As a consequence only a small fraction of previous events can affect the probability 
of future earthquakes so that, after a certain time, Earth crust loses memory of previous seismicity. This aspect is 
perhaps more realistic with respect to the idea, contained in a power law correlation, that events are all correlated 
with each other and also gives rise to important implications for seismic forecasting. The construction of seismic 
catalogs, indeed, dates back to about 50 years, and according to Eq.© one can have good estimates of seismic hazard 
without considering previous seismicity. This is no longer true if one assumes a power law time decorrelation of the 
type © especially for small values of A. We want also to point out that a general state-rate formulation [15j gives rise 
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FIG. 2: (Color online) The intertime distribution obtained using Eq. (|SJ, with two different values of A = 1.2,5 and Ml = 
1.5, 2.5. Continuous and broken curve are the experimental D(At, Ml) with Ml = 1.5 and Ml ~ 2.5 respectively. For A = 1.2 
we set k = 210sec, A = 1.410~ 4 sec _1 , /Lt = 410" 7 , 7 = 1. For A = 5 we set k = 420sec, A = 1.910~ 4 sec~\ fx = 1.510" 6 , 
7 = 0.1. In the inset the magnitude distribution of the experimental catalog (black line) and numerical catalog with A = 1.2 
(red o) and A = 5 (green □). 



to correlations between earthquakes that decay exponentially in time. We finally observe, that taking into account 
only a fraction of previous events in the evaluation of conditional probabilities, the numerical procedure considerably 
speeds up. In the case of long temporal correlation CPU time grows with the number of events as N^, whereas in 
the case of an exponential tail the growth is linear in N e . For this reason, assuming the functional form (JJJ one can 
simulate very large sequences of events. In particular for a different choice of parameters, one can construct synthetic 
catalogs containing the same number of events (N e = 245000 with M > 1.5) of the experimental California Catalog. 
In fig. [3]we compare numerical and experimental distributions D(At, M L ) for three different values of M L . For each 
value of M L , the numerical curve reproduces the experimental data and obviously fulfill Eq.QJ (inset (a) in Fig.{3J)). 
Also the numerical magnitude distribution P(M) is in very good agreement with the experimental one (inset (b) in 
Fig-©)- Finally, evaluation of quantities Qj and Rj for the synthetic catalog leads, as expected, to the same clustering 
behaviour as for experimental data (Fig. lb). After fixing k, we express numerical time unit in seconds and we observe 
that numerical catalog corresponds to a period of about 9.9 * 10 9 sec ~ 30 years. Therefore our model is able to 
construct a synthetic catalog covering about 30 years that contains about the same number of events and displays the 
same statistical organization in magnitude and time of occurrence as real California Catalog. The high efficiency of 
the model in reproducing past scismicity indicates that the model is a good tool for earthquake forecasting. In fact, 
given a seismic history, Eq.J2J together with Eq.s(0][7|) gives the probability to have an earthquake of magnitude M 
at time t inside a considered geographic region. Our approach is different from the Reasenberg- Jones method |l6| . 
which is currently used for evaluation of seismic hazard. This method is based on the generalized Omori law that gives 
for the rate of occurrence of magnitude M aftershocks, P(t, M) — A10 h ( M ~ MM \t — tM + c) - ^, where tM and Mm 
are the time of occurrence and magnitude of the main-shock. The starting set of parameters (b,p, c, A) is estimated 
from previous seismic sequences, and then their value is continuously updated as soon as new data become available. 
However, strong fluctuations in the magnitude distribution observed in Fig.l suggest that the extrapolated b from 
the previous subset may not give the correct value to use for event forecasting. Furthermore, one has an improving 
parameters estimation as the sequence evolves, but at the same time hazard is decreasing. Conversely in our model 
parameters (A, 7) are evaluated on the basis of the entire history of 245,000 events leading to a more precise 
estimation. Nevertheless, due to the stochastic nature of the process, one observes fluctuations of b and p from one 
sequence to the other (Fig. lb). Our model, furthermore, also allows hazard estimation outside the Omori sequence 
and therefore long term forecasting. 

We finally observe that also spatial distributions of seismic events reveal some kinds of scale invariance ^3 ED ■ 
These indicate that also spatial distribution originates from a critical behaviour of the Earth crust suggesting that a 
dynamical scaling hypothesis as in Eq.Q can also work if one appropriately introduces spatial dependencies. In this 
way it would be possible to construct seismic hazard maps. 
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FIG. 3: (Color online) The intertime distribution as a function of A£P c (Ml) obtained using Eq.JTJ (black circle O) and 
compared with the experimental distributions (red □) for three different values of Ml (Ml = 1.5, 2.5, 3.5, from top to bottom). 
We set k — 4.9 * 10 4 sec, A = 6.1 * 10 -5 sec -1 , /i = 2 * 10~ 5 , 7 = 0.1. In inset (a) the scaling behaviour as in Eq.Q and in inset 
(b) the experimental (red) and numerical (black) magnitude distribution. 
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